//
// File: nonlinearHRV.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 23-Apr-2025 16:26:18
//

// Include Files
#include "nonlinearHRV.h"
#include "colon.h"
#include "combineVectorElements.h"
#include "div.h"
#include "freqDomainHRV_data.h"
#include "freqDomainHRV_initialize.h"
#include "freqDomainHRV_types.h"
#include "polyfit.h"
#include "polyval.h"
#include "rt_nonfinite.h"
#include "std.h"
#include "unsafeSxfun.h"
#include "coder_array.h"
#include <cmath>
#include <string.h>

// Function Definitions
//
// nonlinearHRV(IBI,opt) - calculates nonlinear HRV
//
// Inputs:    IBI = inter-beat interval (s) and time locations (s)
//            opt = analysis options from gui
//
// Outputs:   output is a structure containg all HRV.
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
//
//  This file is part of HRVAS
//
//  HRVAS is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  HRVAS is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
// Arguments    : coder::array<double, 2U> &ibi
//                double m
//                double r
//                double n1
//                double n2
//                double breakpoint
//                struct5_T *output
// Return Type  : void
//
void nonlinearHRV(coder::array<double, 2U> &ibi, double m, double r, double n1,
                  double n2, double breakpoint, struct5_T *output)
{
  coder::array<double, 2U> F_n;
  coder::array<double, 2U> c_ibi;
  coder::array<double, 2U> lastrun;
  coder::array<double, 2U> n;
  coder::array<double, 2U> r1;
  coder::array<double, 2U> run;
  coder::array<double, 2U> y;
  coder::array<double, 1U> B;
  coder::array<double, 1U> b_ibi;
  coder::array<double, 1U> b_r;
  coder::array<int, 2U> c_i;
  coder::array<boolean_T, 2U> x;
  double p[2];
  double d;
  double mu;
  double nWin;
  int i;
  int i1;
  int i2;
  int i3;
  unsigned int j;
  int loop_ub_tmp;
  int nj;
  int nx;
  boolean_T exitg1;
  if (!isInitialized_freqDomainHRV) {
    freqDomainHRV_initialize();
  }
  // check inputs
  nj = ibi.size(0) - 1;
  b_ibi.set_size(ibi.size(0));
  for (i = 0; i <= nj; i++) {
    b_ibi[i] = ibi[i + ibi.size(0)] * 1000.0;
  }
  nj = b_ibi.size(0);
  for (i = 0; i < nj; i++) {
    ibi[i + ibi.size(0)] = b_ibi[i];
  }
  // convert ibi to ms
  // assumes ibi units are seconds
  //      if abs(range(ibi(:,2)))<50 %assume ibi units are seconds
  //              ibi(:,2)=ibi(:,2).*1000; %convert ibi to ms
  //      end
  //      if abs(range(diff(ibi(:,1))))>50 %assume time unites are ms
  //          ibi(:,1)=ibi(:,1)./1000; %convert time to s
  //      end
  // function e=sampen(y,M,r);
  //
  // Input Parameters
  //
  //    y:  input signal vector
  //    m:  maximum template length (default m=5)
  //    r:  matching threshold in standard deviations (default r=.2)
  //
  // Output Parameters
  //
  //    e sample entropy estimates for m=0,1,...,M-1
  //
  // Full usage:
  //
  //    [e,se,A,B]=sampen(y,m,r,sflag,cflag,vflag)
  //
  // Input Parameters
  //
  //    sflag    flag to standardize signal(default yes/sflag=1)
  //    vflag    flag to calculate standard errors (default no/vflag=0)
  //
  // Output Parameters
  //
  //    se standard error estimates for m=0,1,...,M-1
  //    A number of matches for m=1,...,M
  //    B number of matches for m=0,...,M-1
  //    (excluding last point in Matlab version)
  //
  //  Ref: http://www.physionet.org/physiotools/sampen/matlab/1.1/
  // normalize ?
  nj = ibi.size(0);
  b_ibi.set_size(ibi.size(0));
  for (i = 0; i < nj; i++) {
    b_ibi[i] = ibi[i + ibi.size(0)];
  }
  nWin = r * coder::b_std(b_ibi);
  // function [e,A,B]=sampenc(y,M,r);
  //
  // Input
  //
  //    y input data
  //    M maximum template length
  //    r matching tolerance
  //
  // Output
  //
  //    e sample entropy estimates for m=0,1,...,M-1
  //    A number of matches for m=1,...,M
  //    B number of matches for m=0,...,M-1 excluding last point
  //
  //  http://www.physionet.org/physiotools/sampen/matlab/1.1/
  i = ibi.size(0);
  nj = ibi.size(0);
  lastrun.set_size(1, ibi.size(0));
  for (i1 = 0; i1 < nj; i1++) {
    lastrun[i1] = 0.0;
  }
  nj = ibi.size(0);
  run.set_size(1, ibi.size(0));
  for (i1 = 0; i1 < nj; i1++) {
    run[i1] = 0.0;
  }
  loop_ub_tmp = static_cast<int>(m);
  b_ibi.set_size(loop_ub_tmp);
  B.set_size(loop_ub_tmp);
  for (i1 = 0; i1 < loop_ub_tmp; i1++) {
    b_ibi[i1] = 0.0;
    B[i1] = 0.0;
  }
  i1 = ibi.size(0) - 2;
  for (int b_i{0}; b_i <= i1; b_i++) {
    nj = (i - b_i) - 2;
    for (loop_ub_tmp = 0; loop_ub_tmp <= nj; loop_ub_tmp++) {
      j = (static_cast<unsigned int>(loop_ub_tmp) + b_i) + 2U;
      if (std::abs(ibi[(static_cast<int>(j) + ibi.size(0)) - 1] -
                   ibi[b_i + ibi.size(0)]) < nWin) {
        d = lastrun[loop_ub_tmp] + 1.0;
        run[loop_ub_tmp] = d;
        i2 = static_cast<int>(std::fmin(m, d));
        for (nx = 0; nx < i2; nx++) {
          b_ibi[nx] = b_ibi[nx] + 1.0;
          if (static_cast<int>(j) < i) {
            B[nx] = B[nx] + 1.0;
          }
        }
      } else {
        run[loop_ub_tmp] = 0.0;
      }
      lastrun[loop_ub_tmp] = run[loop_ub_tmp];
    }
  }
  if (m - 1.0 < 1.0) {
    nj = 0;
  } else {
    nj = static_cast<int>(m - 1.0);
  }
  b_r.set_size(nj + 1);
  b_r[0] = static_cast<double>(ibi.size(0)) *
           (static_cast<double>(ibi.size(0)) - 1.0) / 2.0;
  for (i = 0; i < nj; i++) {
    b_r[i + 1] = B[i];
  }
  B.set_size(b_r.size(0));
  nj = b_r.size(0);
  for (i = 0; i < nj; i++) {
    B[i] = b_r[i];
  }
  // format decimal places
  if (b_ibi.size(0) == B.size(0)) {
    nj = b_ibi.size(0);
    for (i = 0; i < nj; i++) {
      b_ibi[i] = b_ibi[i] / B[i];
    }
  } else {
    rdivide(b_ibi, B);
  }
  nx = b_ibi.size(0);
  for (int k{0}; k < nx; k++) {
    b_ibi[k] = std::log(b_ibi[k]);
  }
  nj = b_ibi.size(0);
  for (i = 0; i < nj; i++) {
    b_ibi[i] = -b_ibi[i] * 1000.0;
  }
  nx = b_ibi.size(0);
  for (int k{0}; k < nx; k++) {
    b_ibi[k] = std::round(b_ibi[k]);
  }
  output->sampen.set_size(b_ibi.size(0));
  nj = b_ibi.size(0);
  for (i = 0; i < nj; i++) {
    output->sampen[i] = b_ibi[i] / 1000.0;
  }
  // DFA(IBI,nn) - calculates alpha from detrended fluctuation analysis
  //
  // Inputs:    y = inter-beat interval signal (s)
  //            n1,n2 = limits of window sizes
  //            breakpoint = value of n that determines where alpha1 ends and
  //            alpha2 begins
  //
  // Outputs:   alpha = slope of log-log plot of integrated y vs window size.
  //
  // Example:   alpha=DFA(y,4,300)
  //
  //  Reference: Heart rate analysis in normal subjects of various age groups
  //  Rajendra Acharya U*, Kannathal N, Ong Wai Sing, Luk Yi Ping and
  //  TjiLeng Chua
  //      [r, c]=size(data);
  //      if r>c; data=data'; end
  if (std::isnan(n1) || std::isnan(n2)) {
    n.set_size(1, 1);
    n[0] = rtNaN;
  } else if (n2 < n1) {
    n.set_size(1, 0);
  } else if ((std::isinf(n1) || std::isinf(n2)) && (n1 == n2)) {
    n.set_size(1, 1);
    n[0] = rtNaN;
  } else if (std::floor(n1) == n1) {
    nj = static_cast<int>(n2 - n1);
    n.set_size(1, nj + 1);
    for (i = 0; i <= nj; i++) {
      n[i] = n1 + static_cast<double>(i);
    }
  } else {
    coder::b_eml_float_colon(n1, n2, n);
  }
  // array of window sizes
  // preallocate memory
  F_n.set_size(1, n.size(1));
  nj = ibi.size(0);
  c_ibi.set_size(1, ibi.size(0));
  for (i = 0; i < nj; i++) {
    c_ibi[i] = ibi[i + ibi.size(0)];
  }
  mu = coder::combineVectorElements(c_ibi) / static_cast<double>(ibi.size(0));
  // mean value
  i = n.size(1);
  for (int b_i{0}; b_i < i; b_i++) {
    double N1;
    d = n[b_i];
    nWin = std::floor(static_cast<double>(ibi.size(0)) / d);
    // number of windows
    N1 = nWin * d;
    // length of data minus rem
    // preallocate memory
    loop_ub_tmp = static_cast<int>(N1);
    lastrun.set_size(1, loop_ub_tmp);
    for (i1 = 0; i1 < loop_ub_tmp; i1++) {
      lastrun[i1] = 0.0;
    }
    // fitcoef=zeros(2,n(i));
    if (N1 < 1.0) {
      nj = 0;
    } else {
      nj = static_cast<int>(N1);
    }
    run.set_size(1, nj);
    for (i1 = 0; i1 < nj; i1++) {
      run[i1] = ibi[i1 + ibi.size(0)] - mu;
    }
    if ((run.size(1) != 0) && (run.size(1) != 1)) {
      i1 = run.size(1);
      for (int k{0}; k <= i1 - 2; k++) {
        run[k + 1] = run[k] + run[k + 1];
      }
    }
    // integrate
    i1 = static_cast<int>(nWin);
    for (nx = 0; nx < i1; nx++) {
      double d1;
      // linear fit coefs
      d = n[b_i];
      nWin = ((static_cast<double>(nx) + 1.0) - 1.0) * d + 1.0;
      d1 = (static_cast<double>(nx) + 1.0) * d;
      if (nWin > d1) {
        i2 = 0;
        i3 = 0;
      } else {
        i2 = static_cast<int>(nWin) - 1;
        i3 = static_cast<int>(d1);
      }
      if (std::isnan(d)) {
        y.set_size(1, 1);
        y[0] = rtNaN;
      } else if (d < 1.0) {
        y.set_size(1, 0);
      } else {
        y.set_size(1, static_cast<int>(d - 1.0) + 1);
        nj = static_cast<int>(d - 1.0);
        for (int k{0}; k <= nj; k++) {
          y[k] = static_cast<double>(k) + 1.0;
        }
      }
      nj = i3 - i2;
      c_ibi.set_size(1, nj);
      for (i3 = 0; i3 < nj; i3++) {
        c_ibi[i3] = run[i2 + i3];
      }
      coder::polyfit(y, c_ibi, p);
      // create linear fit
      d = n[b_i];
      nWin = ((static_cast<double>(nx) + 1.0) - 1.0) * d + 1.0;
      if (nWin > (static_cast<double>(nx) + 1.0) * d) {
        i2 = 1;
      } else {
        i2 = static_cast<int>(nWin);
      }
      if (std::isnan(d)) {
        y.set_size(1, 1);
        y[0] = rtNaN;
      } else if (d < 1.0) {
        y.set_size(1, 0);
      } else {
        y.set_size(1, static_cast<int>(d - 1.0) + 1);
        nj = static_cast<int>(d - 1.0);
        for (i3 = 0; i3 <= nj; i3++) {
          y[i3] = static_cast<double>(i3) + 1.0;
        }
      }
      r1.set_size(1, y.size(1));
      if (y.size(1) != 0) {
        loop_ub_tmp = y.size(1);
        r1.set_size(1, y.size(1));
        for (i3 = 0; i3 < loop_ub_tmp; i3++) {
          r1[i3] = p[0];
        }
        if (y.size(1) == r1.size(1)) {
          nj = y.size(1) - 1;
          r1.set_size(1, y.size(1));
          for (i3 = 0; i3 <= nj; i3++) {
            r1[i3] = y[i3] * r1[i3] + p[1];
          }
        } else {
          binary_expand_op(r1, y, p);
        }
      }
      nj = r1.size(1);
      for (i3 = 0; i3 < nj; i3++) {
        lastrun[(i2 + i3) - 1] = r1[i3];
      }
    }
    //  RMS fluctuation of integraged and detrended series
    if (run.size(1) == lastrun.size(1)) {
      r1.set_size(1, run.size(1));
      nj = run.size(1);
      for (i1 = 0; i1 < nj; i1++) {
        nWin = run[i1] - lastrun[i1];
        r1[i1] = nWin * nWin;
      }
    } else {
      binary_expand_op(r1, run, lastrun);
    }
    F_n[b_i] = std::sqrt(coder::combineVectorElements(r1) / N1);
  }
  // fit all values of n
  x.set_size(1, n.size(1));
  nj = n.size(1);
  for (i = 0; i < nj; i++) {
    x[i] = (n[i] == breakpoint);
  }
  nx = x.size(1);
  nj = 0;
  c_i.set_size(1, x.size(1));
  loop_ub_tmp = 0;
  exitg1 = false;
  while ((!exitg1) && (loop_ub_tmp <= nx - 1)) {
    if (x[loop_ub_tmp]) {
      nj++;
      c_i[nj - 1] = loop_ub_tmp + 1;
      if (nj >= nx) {
        exitg1 = true;
      } else {
        loop_ub_tmp++;
      }
    } else {
      loop_ub_tmp++;
    }
  }
  if (x.size(1) == 1) {
    if (nj == 0) {
      c_i.set_size(1, 0);
    }
  } else {
    if (nj < 1) {
      nj = 0;
    }
    c_i.set_size(c_i.size(0), nj);
  }
  lastrun.set_size(1, c_i.size(1));
  nj = c_i.size(1);
  for (i = 0; i < nj; i++) {
    lastrun[i] = c_i[i];
  }
  // fit short term n=1:bp
  loop_ub_tmp = static_cast<int>(lastrun[0]);
  // fit long term n=bp+1:end
  j = static_cast<unsigned int>(lastrun[0]) + 1U;
  if (j > static_cast<unsigned int>(n.size(1))) {
    i = 0;
    i1 = 0;
  } else {
    i = static_cast<int>(j) - 1;
    i1 = n.size(1);
  }
  if (j > static_cast<unsigned int>(F_n.size(1))) {
    i2 = 0;
    i3 = 0;
  } else {
    i2 = static_cast<int>(j) - 1;
    i3 = F_n.size(1);
  }
  //      lfit=polyval(a,log10(n));
  //      figure; loglog(n,F_n)
  //      hold on; loglog(n,10.^lfit,'r')
  run.set_size(1, n.size(1));
  nj = n.size(1);
  for (int k{0}; k < nj; k++) {
    run[k] = n[k];
  }
  nx = n.size(1);
  for (int k{0}; k < nx; k++) {
    run[k] = std::log10(run[k]);
  }
  y.set_size(1, F_n.size(1));
  nj = F_n.size(1);
  for (int k{0}; k < nj; k++) {
    y[k] = F_n[k];
  }
  nx = F_n.size(1);
  for (int k{0}; k < nx; k++) {
    y[k] = std::log10(y[k]);
  }
  coder::polyfit(run, y, p);
  output->dfa.alpha[0] = std::round(p[0] * 1000.0) / 1000.0;
  output->dfa.alpha[1] = std::round(p[1] * 1000.0) / 1000.0;
  //  total slope
  run.set_size(1, loop_ub_tmp);
  for (int k{0}; k < loop_ub_tmp; k++) {
    run[k] = n[k];
  }
  nj = loop_ub_tmp - 1;
  for (int k{0}; k <= nj; k++) {
    run[k] = std::log10(run[k]);
  }
  y.set_size(1, loop_ub_tmp);
  for (int k{0}; k < loop_ub_tmp; k++) {
    y[k] = F_n[k];
  }
  for (int k{0}; k <= nj; k++) {
    y[k] = std::log10(y[k]);
  }
  coder::polyfit(run, y, p);
  output->dfa.alpha1[0] = std::round(p[0] * 1000.0) / 1000.0;
  output->dfa.alpha1[1] = std::round(p[1] * 1000.0) / 1000.0;
  //  short range scaling exponent
  nj = i1 - i;
  run.set_size(1, nj);
  for (i1 = 0; i1 < nj; i1++) {
    run[i1] = n[i + i1];
  }
  nx = nj - 1;
  for (int k{0}; k <= nx; k++) {
    run[k] = std::log10(run[k]);
  }
  nj = i3 - i2;
  y.set_size(1, nj);
  for (i = 0; i < nj; i++) {
    y[i] = F_n[i2 + i];
  }
  nx = nj - 1;
  for (int k{0}; k <= nx; k++) {
    y[k] = std::log10(y[k]);
  }
  coder::polyfit(run, y, p);
  output->dfa.alpha2[0] = std::round(p[0] * 1000.0) / 1000.0;
  output->dfa.alpha2[1] = std::round(p[1] * 1000.0) / 1000.0;
  //  long range scaling exponent
  output->dfa.F_n.set_size(F_n.size(1));
  nj = F_n.size(1);
  for (i = 0; i < nj; i++) {
    output->dfa.F_n[i] = F_n[i];
  }
  output->dfa.n.set_size(n.size(1));
  nj = n.size(1);
  for (i = 0; i < nj; i++) {
    output->dfa.n[i] = n[i];
  }
}

//
// File trailer for nonlinearHRV.cpp
//
// [EOF]
//
